Journal of Mathematical Biology

Price equation and genotype dynamics are two methods for studying the fixation of one allele by natural selection in a diploid population. There are two strict monotonicity conditions that imply the fixation of one allele. The genotype dynamics is called Haldane monotone if the relative frequency of one allele strictly increases along all solutions of the genotype dynamics, so this allele is fixed. In this paper, we show that the genotype dynamics is Haldane monotone if and only if the right-hand side of the Price equation is always strictly positive. The other strict monotonicity condition requires that the relative frequency of a homozygote strictly increase according to the genotype dynamics. For example, in a model where the genotype dynamics is governed by interactions between individuals, the cost-accepting homozygote is fixed by natural selection if the other genotypes always receive a smaller average gain from all interactions than the cost-accepting homozygote. Both monotonicity conditions require that the interaction is not well-mixed in the population. These two conditions are not equivalent. In addition, we give a non-monotonicity condition, which also implies the fixation of a homozygote. The fixation of a homozygote depends on the phenotypic payoff of the interaction, the genotype-phenotype mapping, and the interaction scheme. In a sexual population, the interaction scheme of siblings depends on the mating system, and so do the conditions of fixation of the cost-accepting homozygote. We present examples showing that if we only change the monogamous mating system, assuming panmixing or mating assortativity, then the condition for the fixation of the cooperator homozygote is b > 2c and b > c, respectively.

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

AO_SCPLOWBSTRACTC_SCPLOWOncolytic viruses provide cancer therapy using replication-competent viruses that selectively infect and lyse tumour cells. Their tumour-specific replication also enables the delivery of targeted, virus-encoded gene products, such as enzymes that activate prodrugs. This dual functionality offers the potential for synergistic effects by combining direct oncolysis with localised drug activation. The interplay between infection, replication, lysis, and gene product delivery remains poorly understood. Here, we introduce a spatially structured, multi-patch model of cancer cells infected by an oncolytic virus engineered to deliver a prodrug-activating enzyme. The spatial system is first represented as a microscopic model and subsequently reduced via spectral dimension reduction techniques. This reduction yields an ordinary differential equation model for a set of coarse-grained variables, which we analyze both without the transgene (OV model) and with the transgene (TOV model). For each case, we compute the basic reproduction number, R0, which determines the conditions for viral spread. Both models exhibit three regimes via transcritical bifurcations: (i) R0 < 0, extinction of both cancer and infected cells; (ii) 0 < R0 [&le;] 1, persistence of cancer cells only; and (iii) R0 > 1, coexistence as a stable node or as a focus. The TOV model, as a difference form the OV model, can undergo periodic oscillations arising from a Hopf-Andronov bifurcation. Notably, oscillation amplitudes can be controlled such that cancer cells largely decrease when drug-induced death is stronger in non-infected cells than in infected ones, enabling effective cancer cells killing while maintaining viral replication and prodrug activation. The qualitative behaviour of the coarse-grained model is shown to be preserved in both the microscopic and spatially explicit models.

We introduce a spectral existence criterion for the evolution of cooperation in the form of the inequality{lambda} maxb > c, where{lambda} max is the leading eigenvalue of an interaction operator encoding population structure, and b and c represent benefit and cost tradeoffs, respectively. Nowaks five rules for the evolution of cooperation correspond to cases in which the cooperation condition reduces to a scalar assortment coefficient. These results follow from the Price equation, which sheds light on a long-standing debate on the role of inclusive fitness and evolutionary dynamics in explaining the evolution of cooperation.

Clavicle fractures often exhibit markedly different clinical outcomes: some patients recover acceptable function despite shortening or displacement, whereas others with apparently similar deformity develop persistent pain, functional loss, or poor healing. To explain this distinction, we propose a minimal nonlinear mechanical model for prognostic analysis of clavicle fractures. The model describes the interaction between fracture-related shortening and compensatory shoulder-girdle posture through a reduced equilibrium equation incorporating stiffness, geometric nonlinearity, and shortening-posture coupling. Within this framework, we analyze equilibrium branches, local stability, and the emergence of critical thresholds. We show that post-fracture destabilization can be interpreted as a fold bifurcation, while more complex parameter dependence gives rise to cusp-type structures and multistability. These bifurcation mechanisms provide a mathematical explanation for sudden deterioration after injury or treatment, as well as for strong inter-individual variability. We further introduce an optimization principle based on a utility functional to guide treatment planning. The analysis predicts that the optimal safe correction should lie strictly below the bifurcation threshold, thereby generating a natural safety margin. Although the model is simplified and has not yet been calibrated against patient data, it nevertheless provides a theoretical framework for understanding why fracture prognosis may deteriorate abruptly near critical mechanical conditions and offers a dynamical-systems interpretation of empirical treatment thresholds used in clinical practice.

The Poisson-Nernst-Planck (PNP) system is an accurate model of electrodiffusion of ionic species. It is commonly used in situations where nanoscale resolution is required, for instance close to ion channels in the membranes of biological cells. The inherent stiffness of the equations has made them challenging to solve and has limited the applicability of the system. In particular, the time step required for stable solutions has typically needed to be very short (nanoseconds), which makes simulations on the time scale of an action potential (milliseconds) difficult. Recently, it has been observed that avoiding operator splitting and instead solving the concentration equations and the electrostatic equation in a coupled manner relaxes the time-step limitation considerably. However, no theoretical explanation of this observation has been provided. Here, we aim to explain why the coupled scheme allows much larger time steps. We illustrate the mechanism by considering special cases that define necessary, but not sufficient, conditions for stability. We also show that these conditions remain relevant for the fully coupled PNP model in 3D.

1Although resources are typically distributed continuously in space, species distributions often organize into discrete clusters. In his seminal paper [36], Turing demonstrated that such clusters can spontaneously arise in population densities, even when populations evolve in environments with continuously varying conditions. This phenomenon is known as Turing instability. In this work, we focus on two models grounded in population dynamics: a one-dimensional model based on the nonlocal Fisher-KPP equation, and a two-dimensional model involving an environmental gradient. We show that phenotypic clusters (sometimes referred to as "species") emerge in these models. We prove that they do not emerge because of Turing instability, but because of stochasticity, and that they disappear when stochasticity is reduced. First, for both models, we start our simulations with initial populations uniformly distributed in the state space. We show that phenotypic clusters quickly emerge and that the distances between them depend on the population size, that is, on the degree of stochasticity. Next, we start from already clearly defined phenotypic clusters. We identify three regimes in the connection between population size, the initial distances between clusters, and the distances between clusters at equilibrium. Last, on the two-dimensional model, we relax the hypothesis of complete clonality by varying the effective recombination rate, explore its effect on phenotypic clustering, and show that phenotypic clustering decays drastically with slight recombination.

Experimental studies show that tumor cells adopt migratory or proliferative phenotypes depending on the local extracellular matrix (ECM). In this work, we propose a minimal go-or-grow invasion model, comprising two specialist cell phenotypes: proliferating and migratory, with phenotypic switching and cell migration depending on local ECM density. Numerical simulations of this model reveal that the spatial arrangement of proliferative and migratory cells depends on the choice of phenotypic switching function. We then ask whether this specialist cell-phenotype model can be reduced to a generalist cell-phenotype model. We derive a relationship between the reduced model and go-or-grow model in the fast phenotypic switching regime. We observe that the reduced model captures the dynamics of the original model, for a range of realistic phenotypic switching functions. We analytically derive the minimum traveling wave speed of the reduced model in a homogeneous ECM bed. Moreover, using linear stability analysis on the go-or-grow model, we recover the same wave speed expression. In addition, we numerically explore how the key parameters influence the traveling wave speed profile. Our analysis indicated the counter-intuitive result that the wave speed is independent of the matrix degradation rate, and our simulations show that, at most, the speed is weakly dependent on this parameter.

This work focuses on the global and partial identification problem for fractional differential equations. We provide a general numerical procedure based on global and local optimization algorithms with two refinements for biological systems that ensure solution positivity and homogeneous parameter units. The method is applied to a new fractional model of Dengue outbreak called the Fractional Homogeneous Nishiura (FHN) model, calibrated using data of newly infected people in Cape Verde. We show that our identification method yields a better fit between data and model solutions than previous approaches and that our FHN model captures the dynamics of Dengue more closely than existing systems.

Biological tissues grow at rates that depend on the geometry of the supporting tissue substrate. In this study, we present a novel discrete mathematical model for simulating biological tissue growth in a range of geometries. The discrete model is deterministic and tracks the evolution of the tissue interface by representing it as a chain of individual cells that interact mechanically and simultaneously generate new tissue material. To describe the collective behaviour of cells, we derive a continuum limit description of the discrete model leading to a reaction-diffusion partial differential equation governing the evolution of cell density along the evolving interface. In the continuum limit, the mechanical properties of discrete cells are directly linked to their collective diffusivity, and spatial constraints introduce curvature dependence that is not explicitly incorporated in the discrete model. Numerical simulations of both the discrete and continuum models reproduce the smoothing behaviour observed experimentally with minimal discrepancies between the models. The discrete model offers further individual-level details, including cell trajectory data, for any restoring force law and initial geometry. Where applicable, we discuss how the discrete model and its continuum description can be used to interpret existing experimental observations.

Cell-cell adhesion is a key regulator of cancer invasion. In this work, we extend a pre-existing individual based cancer invasion model by introducing a stochastic representation of N-cadherin-mediated adhesion, where the lifetime of a cell-cell bond depends on the pulling force acting on the bond. Using experimental data, we derive expressions for the mean and standard deviation of N-cadherin bond lifetimes and fit them to Gamma distributions, enabling their treatment as force-dependent random variables. These distributions are then used to modify the diffusion coefficient of mesenchymal cancer cells. The model predicts reduced random motility with increasing adhesion and incorporates a dynamic transition between catch- and slip-bond behaviour. Along with this model for cell motility, we propose a preliminary physical framework, that can be used to model pattern formation as a result of the new adhesion mechanic.

The influence of arbitrary randomness in cell division times on the variability of protein copy numbers within a lineage ensemble has been recently studied, going beyond the contributions of noisy gene expression and partitioning error. However, variability of protein concentrations need separate study, since cell size growth between cell divisions dilute protein concentrations at the same rate as size growth, which also determines mean division times. Here for a model of bursty protein production, we present exact moments (of all orders) of protein concentrations in the cyclo-stationary state, comparing: (i) population and lineage cell ensembles, and (ii) statistics at different cell ages. Two interesting results emerge. While the variance of protein concentration changes with the degree of division time heterogeneity at any cell age, the age-averaged variance is independent of it within lineage ensemble but stays dependent within population ensemble. The skewness within population ensemble is higher in younger cells than within lineage ensemble, and this behavior reverses at older ages. Such a feature vanishes for the age-averaged distribution, with population based skewness always dominating over that of lineage. We also show that mother-daughter correlations in generation times, do not add any significant difference to the results.

Spatio-temporal interactions shape microbial community dynamics. Metabolism, through competition and cross-feeding, is a foundational mechanism of these interactions. Flux balance analysis enables efficient simulation of steady-state metabolism. Integrating these simulations through time, using dynamic flux balance analysis, provides temporal predictions of growth and metabolism. Incorporating spatial context, through partial differential equations, enables spatio-temporal simulation of microbial communities. In this chapter, we step through this sequential process, moving from steady-state, to temporal, to spatio-temporal simulation of microbial community metabolism. We provide an illustrative example using the modeling software COMETS (Computation of Microbial Ecosystems in Time and Space) to simulate interacting bacterial colonies of Bifidobacterium longum subsp. infantis and Anaerobutyricum hallii (previously Eubacterium hallii). Within this simulation, both competition and cross-feeding influenced the production of butyrate leading to an intermediate optimal interaction distance for metabolite production. We outline each step and provide open-source code such that this simulation can serve as a template for future spatio-temporal simulations of microbial community metabolism.

Oscillatory behaviour is important in multiple biological contexts. However, the inherent nonlinearity and high dimensionality of mathematical models in biology makes proving the existence and the localization of limit cycle oscillations challenging. Here, we provided an elementary proof for the existence and a method for rigorously localizing the oscillatory solutions in a class of benchmark biomolecular oscillators. To construct the proof, we used a geometric approach based on Brouwers Fixed Point theorem. We constructed a toroidal-like manifold within a positively invariant set by removing the hypervolume containing the fixed point and the trajectories converging to it along its stable manifold. We showed that the vector field describing the system dynamics maps a cross section of the toroidal-like manifold onto itself. The existence of a limit cycle solution in this manifold was guaranteed by Brouwers Fixed Point theorem. For different sets of initial conditions in these cross-sections, we used an interval-based Reachability Analysis to localize the oscillatory behaviour that complements the Brouwers Fixed Point theorem approach. These results add a simple and elegant approach to demonstrating the existence of limit cycles in biomolecular systems as well as a method for rigorous localization of the region of existence.

Minimum viable populations (MVPs) are population levels large enough to surmount risk from demographic, environmental, and genetic stochasticity. MVPs are estimated by biologists to guide conservation practices. However, MVPs are generally estimated for a target population without regard for how they interact with intra- and inter-species population dynamics in the broader ecological community. Thus, how and why population dynamics interact with MVPs imposed by conservation biologists remain unclear. When MVPs are imposed on a continuous population model, traditional analyses fail to capture the range of possible outcomes those MVPs create. Here, we describe viability space decomposition (VSD) as a mathematical tool to systematically analyze the potential crossing of MVPs during population dynamics. We demonstrate that different extinction and survival outcomes can be recovered from a model with imposed MVPs using three VSD concepts in junction with a traditional phase portrait: mortality manifolds which separate conditions that lead to different existential outcomes, ordering manifolds which determine the order of extinction events for multiple populations, and collapse manifolds which determine the survival or extinction of one species given the loss of another. We employ these methods with a standard consumer-resource model, and the methods can be scaled to systems with more species. VSD is a useful tool for conservation biologists and community ecologists concerned with boundary crossing problems in any dynamical system.

Darwinian fitness is equated here with invasion fitness and defined as the quantity determining the fate--certain extinction or possible spread--of a single mutant type. We derive it, together with its phenotypic derivative, for evolution in group-structured populations under limited genetic mixing, where the demography of the focal species and its environment is modeled as a discrete-time stochastic process. Reproduction, physiological development, dispersal, and survival are influenced by interactions within and between groups and by environmental fluctuations within and across generations. Using multitype branching processes in random environments, we show that invasion fitness is predicted by a stochastic growth rate that can be represented biologically in two meaningful genealogical ways. First, as the long-term geometric mean of the expected per-capita number of mutant copies produced per time step by a representative member of the mutant lineage. Second, as the the long-term geometric mean of the expected reproductive-value-weighted per-capita number of mutant copies produced by such an individual. This latter representation is useful for computing the phenotypic directional derivative of invasion fitness. Moreover, this derivative can be written as an actor-centered inclusive-fitness effect derived from properties of the resident population process. This effect depends on class-specific fitness differentials, relatedness, reproductive values, and class frequencies. However, unless generation- and class-specific fitness defines a stochastic matrix, the derivative does not separate stochastic reproductive values from relatedness and class frequencies, and must be evaluated by simulations. In summary, we formalize invasion fitness biologically quite generally and show how Hamiltons marginal rule is deduced from it.

The basic reproduction number of an infectious disease is known to depend on the structure of contacts between individuals in a population. This relationship has been explored mathematically through two well-known models: one which depends on a matrix of contact rates between different demographic groups, and another which depends on the variability of contact rates over the population. Here we introduce a model that combines and generalises these two approaches. We derive a formula for the basic reproduction number and validate it through comparisons to simulated outbreaks. Applying this method to contact survey data collected in Belgium between 2020 and 2022, we find that our model produces higher estimates of the basic reproduction number and larger relative changes over periods when social contact behaviour was changing during the COVID-19 pandemic. Our analysis suggests some practical considerations when using contact data in models of infectious disease transmission.

We introduce a mechanistic, nonlocal tumour-growth model designed specifically to capture explosive dynamics that are not adequately explained by standard logistic reaction-diffusion descriptions. The motivation is empirical: the universal scaling law reported in [1] provides compelling cross-sectional evidence of superlinear tumour activity versus tumour burden, but as a phenomenological relationship it does not by itself supply a dynamical mechanism, nor does it rigorously describe how explosive growth emerges, how fast it develops, or how spatial interactions and tissue boundaries influence it. Our model addresses this gap by incorporating nonlocal proliferative feedback--cells respond to a spatially aggregated neighbourhood signal--and a singular, Kawarada-type acceleration that produces "quenching": tumour density stays bounded while the proliferative drive becomes unbounded as the aggregated signal approaches a critical threshold. This offers a concrete mechanistic route to explosive escalation consistent with physical boundedness. We analyse the model under no-flux (Neumann) boundary conditions, appropriate for reflecting tissue interfaces. In the spatially homogeneous setting we prove finite-time onset of the explosive regime and obtain explicit rates for how rapidly it is approached. For spatially heterogeneous perturbations we derive a transparent spectral stability theory showing how the interaction kernel selects spatial scales and how the singular acceleration tightens stability margins as the explosive threshold is approached. These results provide interpretable links between nonlocal interaction structure, boundary effects, and the emergence of rapid growth. Finally, to connect mechanism to data in the spirit of [1], we embed the model in a Bayesian inference framework that treats the interaction kernel and the acceleration strength as unknown and learned from tumour-growth observations. This enables uncertainty-aware estimation of explosive onset times, escalation rates, and stability margins, while positioning the scaling law of [1] as an observable signature that our mechanistic model can explain and quantify rather than merely fit.

Cancer-associated fibroblasts (CAFs) are a central component of the tumor microenvironment that facilitate a supportive niche for cancer progression and metastasis. Experimental evidence suggests that CAFs can facilitate estrogen-independent tumor growth, thereby reducing the efficacy of anti-hormonal therapies. Understanding and quantifying the complex interactions between tumor cells, hormonal signalling, and the microenvironment are crucial for designing more effective and individualized treatment strategies. We propose a mathematical framework to explore the influence of CAFs on ER+ breast cancer progression and to evaluate strategies to mitigate their impact. We develop a system of nonlinear ordinary differential equations that substantiates the experimental observations by providing a mechanistic basis for the role of CAFs in regulating estrogen-independent growth dynamics. We then employ optimal control theory to evaluate distinct therapeutic approaches involving monotherapy or combinations of: (i) inhibition of tumor-to-CAF signaling, (ii) inhibition of CAF-to-tumor proliferative signaling, and (iii) endocrine therapy. Taken together, our results demonstrate that CAF-targeted strategies can enhance treatment efficacy across various estrogen dosing regimens. Our study provides new insights into the potential of CAF as a therapeutic target that could help to improve existing approaches for endocrine therapies.

Variability in gene expression among single cells and growing cell populations can arise from the stochastic nature of protein synthesis, which often occurs in random bursts. This study investigates the variability in the expression of a growth-sustaining protein, whose concentration is regulated by a negative feedback loop due to cell growth-induced dilution. We model the distribution of protein concentration using a Chapman-Kolmogorov equation for single cells and a population balance equation for growing cell populations. For single cells, we derive an explicit solution for the protein concentration distribution in state space and represent it as a Bessel function in Laplace space. For growing populations, we find that the distribution satisfies a Heun differential equation with singular boundary conditions. By addressing the central connection problem for the Heun equation, we quantify the population-level protein distribution and determine the Mathusian parameter, which characterizes population growth. This work provides a comprehensive analytical framework for understanding how stochastic protein synthesis impacts gene expression variability and population dynamics.